!\Description:
!  Generate a random initial residual vector for the Arnoldi process.
!  Force the residual vector to be in the range of the operator OP.
!\Arguments
!  IDO     Integer.  (INPUT/OUTPUT)
!          Reverse communication flag.  IDO must be zero on the first
!          call to cgetv0.
!          -------------------------------------------------------------
!          IDO =  0: first call to the reverse communication interface
!          IDO = -1: compute  Y = OP * X  where
!                    IPNTR(1) is the pointer into WORKD for X,
!                    IPNTR(2) is the pointer into WORKD for Y.
!                    This is for the initialization phase to force the
!                    starting vector into the range of OP.
!          IDO =  2: compute  Y = B * X  where
!                    IPNTR(1) is the pointer into WORKD for X,
!                    IPNTR(2) is the pointer into WORKD for Y.
!          IDO = 99: done
!          -------------------------------------------------------------
!
!  BMAT    Character*1.  (INPUT)
!          BMAT specifies the type of the matrix B in the (generalized)
!          eigenvalue problem A*x = lambda*B*x.
!          B = 'I' -> standard eigenvalue problem A*x = lambda*x
!          B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
!
!  ITRY    Integer.  (INPUT)
!          ITRY counts the number of times that cgetv0 is called.
!          It should be set to 1 on the initial call to cgetv0.
!
!  INITV   Logical variable.  (INPUT)
!          .TRUE.  => the initial residual vector is given in RESID.
!          .FALSE. => generate a random initial residual vector.
!
!  N       Integer.  (INPUT)
!          Dimension of the problem.
!
!  J       Integer.  (INPUT)
!          Index of the residual vector to be generated, with respect to
!          the Arnoldi process.  J > 1 in case of a "restart".
!
!  V       Complex N by J array.  (INPUT)
!          The first J-1 columns of V contain the current Arnoldi basis
!          if this is a "restart".
!
!  LDV     Integer.  (INPUT)
!          Leading dimension of V exactly as declared in the calling
!          program.
!
!  RESID   Complex array of length N.  (INPUT/OUTPUT)
!          Initial residual vector to be generated.  If RESID is
!          provided, force RESID into the range of the operator OP.
!
!  RNORM   Real scalar.  (OUTPUT)
!          B-norm of the generated residual.
!
!  IPNTR   Integer array of length 3.  (OUTPUT)
!
!  WORKD   Complex work array of length 2*N.  (REVERSE COMMUNICATION).
!          On exit, WORK(1:N) = B*RESID to be used in SSAITR.
!
!  IERR    Integer.  (OUTPUT)
!          =  0: Normal exit.
!          = -1: Cannot generate a nontrivial restarted residual vector
!                in the range of the operator OP.
!-----------------------------------------------------------------------
subroutine cgetv0(ido, initv, n, j, v, resid, rnorm, workd, ierr)
    implicit none
    logical    initv
    integer    ido, ierr, j, n
    real       rnorm
    complex    resid(n), v(n, j), workd(2 * n)
    ! Parameters
    complex, parameter :: one = (1.0, 0.0), zero = (0.0, 0.0)
    complex, external :: cdotc
    real, external :: scnrm2, slapy2
    ! Local Scalars & Arrays
    logical    first, inits, orth
    integer    idist, iseed(4), iter
    real       rnorm0
    save       first, iseed, inits, iter, orth, rnorm0
    ! Data Statements
    data       inits /.true./

    ! Initialize the seed of the LAPACK random number generator
    if (inits) then
        iseed(1) = 1
        iseed(2) = 3
        iseed(3) = 5
        iseed(4) = 7
        inits = .false.
    end if

    if (ido == 0) then
        ierr = 0
        iter = 0
        first = .FALSE.
        orth = .FALSE.
        !-----------------------------------------------------%
        ! Possibly generate a random starting vector in RESID |
        ! Use a LAPACK random number generator used by the    |
        ! matrix generation routines.                         |
        !    idist = 1: uniform (0,1)  distribution;          |
        !    idist = 2: uniform (-1,1) distribution;          |
        !    idist = 3: normal  (0,1)  distribution;          |
        !-----------------------------------------------------%
        if (.not. initv) then
            idist = 2
            call clarnv(idist, iseed, n, resid)
        end if
    end if
    ! Back from computing B*(initial-vector)
    if (first) go to 20
    ! Back from computing B*(orthogonalized-vector)
    if (orth)  go to 40
    ! Starting vector is now in the range of OP; r = OP*r;
    ! Compute B-norm of starting vector.
    first = .TRUE.
    call ccopy(n, resid, 1, workd, 1)
20  continue

    first = .FALSE.
    rnorm0 = scnrm2(n, resid, 1)
    rnorm = rnorm0
    ! Exit if this is the very first Arnoldi step
    if (j == 1) go to 50

    !----------------------------------------------------------------
    ! Otherwise need to B-orthogonalize the starting vector against |
    ! the current Arnoldi basis using Gram-Schmidt with iter. ref.  |
    ! This is the case where an invariant subspace is encountered   |
    ! in the middle of the Arnoldi factorization.                   |
    !                                                               |
    !       s = V^{T}*B*r;   r = r - V*s;                           |
    !                                                               |
    ! Stopping criteria used for iter. ref. is discussed in         |
    ! Parlett's book, page 107 and in Gragg & Reichel TOMS paper.   |
    !---------------------------------------------------------------%
    orth = .TRUE.
30  continue

    call cgemv('C', n, j - 1,  one, v, n, workd, 1, zero, workd(n + 1), 1)
    call cgemv('N', n, j - 1, -one, v, n, workd(n + 1), 1, one, resid, 1)
    ! Compute the B-norm of the orthogonalized starting vector
    call ccopy(n, resid, 1, workd, 1)
40  continue

    rnorm = scnrm2(n, resid, 1)
    ! Check for further orthogonalization.
    if (rnorm > 0.717 * rnorm0) go to 50

    iter = iter + 1
    if (iter <= 1) then
        ! Perform iterative refinement step
        rnorm0 = rnorm
        go to 30
    else
        ! Iterative refinement step "failed"
        resid(1:n) = 0.0
        rnorm = 0.0
        ierr = -1
    end if
50  continue
    ido = 99
end
